function H = CalculateHomography3x3 (im1Coords, im2Coords)

% solution least squares
% A = [x, y, 1, x, y, 1]
% x = [a,b,c,d,e,f]
% y = [x'+y']
% A'y=A'Ax
% need to solve x, knowing A and y
vecOnes = ones(1,length(im1Coords(1,:)));
vecZeros = zeros(1,length(im1Coords(1,:)));
A = [im1Coords(1,:), vecZeros; im1Coords(2,:), vecZeros; vecOnes, vecZeros; ...
    vecZeros, im1Coords(1,:); vecZeros, im1Coords(2,:); vecZeros, vecOnes]';
y = [im2Coords(1,:), im2Coords(2,:)]';
AA = A'*A;
yy = A'*y;
x = AA^(-1)*yy

H = [x(1), x(2), x(3); x(4), x(5), x(6); 0, 0, 1];
